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Abstract 

Linearized numerical stability bounds for solving the nonlinear time-dependent Schrodinger equa- 
■^r • tion (NLSE) using explicit finite-differencing are shown. The bounds are computed for the fourth-order 

Runge-Kutta scheme in time and both second-order and fourth-order central differencing in space. Re- 
sults are given for Dirichlet, modulus-squared Dirichlet, Laplacian-zero, and periodic boundary conditions 
y^ • for one, two, and three dimensions. Our approach is to use standard Runge-Kutta linear stability theory, 

' treating the nonlinearity of the NLSE as a constant. The required bounds on the eigenvalues of the 

scheme matrices are found analytically when possible, and otherwise estimated using the Gershgorin 
CO ■ circle theorem. 

>: 
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'^. ■ 1 Introduction 

•>: 

^^ . The nonlinear Schrodinger equation (NLSE) is used to model a wide variety of physical systems since it 

describes, to least nonlinear order, modulated wave propagation [1]. The general form of the NLSE can be 
written as 

i^ + aV2^-y(r)* + s|*p* = 0, (1) 

X. ... 

^ I where 4* G (D is the value of the wavefunction, V^ is the Laplacian operator, and where a > and s are 

parameters defined by the system being modeled. V{r) is an external potential term, which when included, 
makes Eq. ([1]) known as the Gross-Pitaevskii equation ''2[ . 

Often, efficient and easy-to-use numerical methods are employed to simulate the NLSE. One such method 
is the method of lines where the time-stepping and spatial differencing are treated independently. This trans- 
forms the partial differential equation (PDE) into a large number of coupled ordinary differential equations 
(ODEs). These ODEs can then be solved using a variety of numerical schemes, one of the most common 
being the fourth order Runge-Kutta (RK4) scheme [3]. Using the RK4 scheme with the NLSE produces a 
fully explicit scheme where each grid point at time t is only a function of values at time t — k where k is the 
time-step. This simplifies computational implementations because no matrices are needed to be formed and 
stored, and no linear systems are needed to be solved (which in the nonlinear case also require a nonlinear 
iterative process). 
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The only drawback to using explicit finite-difference schemes (such as the RK4) for simulating PDEs 
is that they are conditionally stable. This means that there is an upper bound on the allowed size of the 
time-step which is dependent on the spatial-step size. If the time-step is larger than this bound, the scheme is 
unstable and diverges [1] . Although rough estimates of the stability bound can be found through an inefficient 
educated guess-and-check, for higher dimensional scenarios, as well as long and/or large simulations, a more 
refined and predictable stability bound is essential for efficient simulations. 

In this paper, we formulate linearized stability bounds for simulating the NLSE with the RK4 scheme. 
The stability bounds depend on the specific spatial differencing scheme being used, as well as on the bound- 
ary conditions. We formulate the bounds for both second-order and fourth-order spatial differencing with 
a variety of boundary conditions (Dirichlet, modulus-squared Dirichlet (MSD), Laplacian-zero (LO), and 
periodic). Each analysis is done for one, two and three dimensions. 

The paper is organized as follows. In Section[2]we review the basic RK4 stability properties and apply the 
results to the NLSE to formulate general stability bounds. Our basic procedure in finding the stability bounds 
and the linear algebra theorems that we utilize are also discussed. In Section [3] we summarize the forms 
of the boundary conditions we consider. Our main analysis begins in Section 3] with the one-dimensional 
NLSE. Linearized stability bounds are found for each scheme and boundary condition combination. In 
Sections [S] and [HI we use the same procedures to formulate the bounds for the two- and three-dimensional 
NLSE respectively. In Section [3 we conclude and summarize all the results from Sections HI [5l and [6] into a 
concise reference. 

2 Stability Theory 

2.1 General Runge-Kutta Scheme Stability 

Given an initial value problem of a set of linear first-order ODEs (in our case, a method-of-lines explicit 
PDE finite-difference scheme), one can formulate the matrix notation 

^-A% (2) 

where A contains the coefficients of the right-hand-sides of the ODEs. We now define 

p = kA, (3) 

where k is the time-step size and A contains the eigenvalues of A. In our case, the eigenvalues of A will 
have the spatial-step size (denoted h) included in them, as well as any parameters of the NLSE. As shown 
in Ref. [F, for the fourth-order Runge-Kutta scheme, if a vector R{p) is defined whose elements are the 
polynomials 

then the stability of the RK4 scheme is guaranteed if 

P(p)ll=o < 1, (5) 

where || |]oo denotes the infinity norm defined as ||a;||oo = max{|a;o|, \xi\, ..., |a;Ar_i|}. Inserting Eq. ^ into 
Eq. (gl) yields 

|i?(A)r = l + ^k«|A|«~lk«|Ar-f(^-k4|A|^ + 24)iiRe(A) (6) 

+ (k^|A|^ + 24) ^ (Re(A))2 + (k^|Ap + 4) ^ (Re(A))3 + ^ (Re(A))^ 
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Figure 1: (Color online) Left: Stability regions for Runge-Kutta schemes. Schemes from first-order to 
fourth-order arc shown from center outwards. Right: Magnified view of the same plot near the point where 
Re(A) = 0. 



In Fig. [T] we show the stability region for the RK4 scheme given by Eq. ^ as well as that for lower-order 
Runge-Kutta schemes (whose R{p) is defined by progressively truncated versions of Eq. (j4])) [5]. As we shall 
show, the eigenvalues of the A matrix are all purely imaginary (or nearly so) in the case of the nonlinear 
(and linear) Schrodinger equation. Thus, it can be seen from Fig. [1] that the third-order Runge-Kutta is 
the lowest-order RK scheme that is conditionally stable for the Schrodinger equations (however, as shown 
in Ref. [6^, this is not the case if the real and imaginary parts of the NLSE are computed in a staggered 
time grid, or, as in Ref. i?;, an artificial dissipative term is added to the NLSE), with the RK4 yielding a 
significantly larger bound on k. This is in contrast to similar PDEs such as the heat equation, whose A 
matrix eigenvalues are typically all real- valued, in which case even forward differencing (RKl) is conditionally 
stable. 

If, as in our case, Re(A) = 0, Eq. ^ simplifies greatly and becomes 
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in which case, Eq. ([5]) leads to the simple stability bound 



k< 



V8 

A||oc 



(7) 



(8) 



2.2 Application to the NLSE 

Applying the above stability theory to the NLSE has the obvious problem that the analysis is purely linear, 
while the NLSE has one (or more) nonlinear terms. A full nonlinear stability analysis is beyond the scope of 
this paper, so instead we linearize the problem by treating the nonlinearity (|^P) as a constant value U (the 
external potential term is usually a constant independent of ^ at each grid point, and so does not need any 
special treatment). This has been done previously for the one-dimensional coupled NLSE for fourth-order 
differencing (in the exclusive case where s < 0) in Ref. ^. Since the value of |^p changes over time during 
the simulation, the linearized stability bound will also change over time. This change in many cases is 
expected to be small (which we have confirmed in numerical simulations, not reported here) and therefore 



can be ignored, i.e. one may compute the bound using the initial condition of ^ (and V{r)) and just leave 
a few percent leeway to cover any changes. This is especially true in the repulsive case (s < 0) where most 
situations have a constant-density background (or maximum background) and the dynamics do not cause 
the maximum background value to change significantly (for example, when simulated coherent structures, 
most of the dynamics are translations of the initial condition with little change in structure). In attractive 
cases (s > 0), blow-up can occur which can alter the stability bound greatly, causing the simulation to crash 
(although in such a case the wavefunction is exploding towards infinity, which most finite-difference schemes 
cannot handle anyways). Many times, simulations of a steady-state or near-steady-state in the modulus- 
squared with a constant potential are performed. In such situations, the linearized stability bounds will be 
(nearly) exact. 

It is also useful to formulate stability bounds for the linear Schrodinger equation (LSE) (where s = 
and V{r) =0). In addition to providing bounds for the LSE, as will be discussed below, the results can also 
be used as practical estimates of the stability bounds for the NLSE (the discrepancy can often be solved by 
lowering the bound by a few percent). 

2.3 Stability Analysis Procedure 

In order to simplify the analysis, we first rewrite Eq. ([2]) as 

— ~ A'i' = — A^ 

at h^ ' 

where h is the step-size of the spatial finite-difference scheme being used. Then, assuming all eigenvalues of 
A are real- valued, the stability condition of Eq. ^ becomes 

PaIIoo a 

In order to be able to use the stability bound of Eq. ©, we must first confirm that all eigenvalues 
of A are purely real (or nearly so) for each scheme/boundary condition combination. In cases where the 
eigenvalues arc not able to be easily computed analytically, we show that the A matrix's eigenvalues are a 
set of boundary values with the remaining eigenvalues being those of a symmetric matrix denoted A' . Then 
by Thm. [l] it is known that all the eigenvalues of A are real. 

Theorem 1 (Ref. [9^ ) . The eigenvalues of a real symmetric matrix are real. 

Once it has been established that Eq. ([9]) can be used, in order to get an upper-bound on k, we require an 
upper-bound on the maximum absolute eigenvalue of A. Due to the sparsity and diagonal dominance of A, 
a good estimate of the upper-bound can be found using the Gershgorin circle theorem (Def. [1] and Thm. [2]) . 

Definition 1 (Ref. jlOl). Let A be a square complex matrix. Around every element an on the diagonal of the 
matrix, a circle with radius equal to the sum of the norms of the other elements in the same row (X^iatj I'^ijl) 
is known as a Gershgorin disc. 

Theorem 2 (Ref. [10 ). Every eigenvalue of a square complex matrix A lies in one of its Gershgorin discs. 

Since every eigenvalue must be contained in a Gershgorin disk, by finding the maximum absolute value 
of the limits of the disks will yield an upper-bound on the maximum modulus of the eigenvalues of A. 

In the one-dimensional LSE case with no external potential and periodic boundary conditions, the A 
matrix becomes circulant as defined by Def. [21 In this case, the eigenvalues can be computed analytically 
by Thm. [31 The upper-bound is then taken by finding the limit of the maximum eigenvalue as the size of 
the matrix goes to infinity. 

Definition 2 (Ref. [11]). A circulant matrix is a square N x N matrix C that can be fully specified by one 
vector, c — {co,Ci, ...,cn-i}, which appears as the first column of C . The remaining columns of C are each 
cyclic permutations of the vector with the offset equal to the column index. 



Theorem 3 (Ref. [11'). The eigenvalues of a circulant matrix are given by 

\j = Co + CN-lUJj + CN~2^] + ■■■ + CiUjf^^, j = 0, ..., A^ 



2Tr i j 

N 



where 

ujj ~ exp 

3 Boundary Conditions 

Since boundary conditions of the spatial differencing in a PDE like the NLSE have the potential to alter 
the stability of a scheme, it is necessary to have stability results for each specific boundary condition one 
would like to use. In this paper we limit ourselves to four boundary conditions which we feel are a good 
combination of simplicity and usefulness. These boundary conditions are: periodic, Dirichlet, Laplacian-zero, 
and Modulus-Squared-Dirichlet. As notation, we use the subscript b to represent any boundary point, and 
b — 1 to represent the grid position one point inward from the boundary in the normal direction. 

For use with the stability analysis, it is desirable to formulate each boundary condition in terms of the 
temporal derivative in the form 

'"i?b*b, (10) 



'dt 
and in terms of the spatial Laplacian in the form 



h^' 



\7Hb = ^Db^b, (11) 

where Bb and Db are assumed to be real- valued constants (possibly differing per boundary point) and defined 
based on the specific boundary condition being used. For periodic boundary conditions (or linear one-sided 
conditions not discussed here), these forms are not applicable. Writing the boundary conditions in the forms 
of Eq. (ITUl) and Eq. pT|) allows them to be expressed in the A matrix as a single real- valued entry (Bb), and 
in the case of the form of the fourth-order differencing chosen here (see Sec. l4.2p . the near-boundary interior 
points will contain Db in their formulation. 

3.1 Periodic 

For periodic boundary conditions, any element of the scheme that is too small or too large in index (i.e. they 
are 'off the grid') are simply replaced by the grid points on the opposite side of the grid. In the case of the 
NLSE, periodic boundary conditions can be problematic especially in background-density situations due to 
the unpredictable phase jump from one side of the grid to the other. 

3.2 Dirichlet 

Dirichlet boundary conditions are defined as 

where C is a constant. In terms of the temporal derivative of the NLSE, this condition is 



5* 

'dt 



0, 



in which case Bb = in Eq. ()T0|) . When inserted into the NLSE, this condition in terms of the Laplacian is 
given by 

VHb = --{s\^b\^-Vb)^b, 
a 

and therefore Db = -/iVa(s|*bP - H) in Eq. i^^. 



3.3 Modulus-Squared Dirichlet 



In some situations Dirichlet boundary condition can fail. Such failure typically occurs in simulations with a 
constant-density background, i.e. a constant value of j'i'p at the boundaries. A standard Dirichlet condition 
will not work in such cases because it does not take into account the phase rotation of '5. Instead, one would 
like to have the modulus-squared of the wavefunction to be constant at the boundaries, i.e. 



1*6 



= C, 



where C is a constant. We have recently formulated a method for such a boundary condition (which is almost 
as easy to implement as Dirichlet) called the modulus-squared Dirichlet boundary condition [T^] . The MSD 
boundary condition is given in terms of the temporal derivative of the NLSE as 



* 



t,b ■ 



zim 



* 



t,b-l 



* 



6-1 



*6 



(12) 



where d'^b-i/dt is computed by the interior scheme first, and then used to compute the boundary values. 



Using the MSD boundary condition gives Bb = 

constant independent of *. As shown in Ref. [H 

condition, Eq. (rT2|) can be viewed as 

a* 
'dt 



{h^/a)lm 



dt \b- 



which is nonlinear, and not a 



'1 *b-i 
due to the underlying assumptions of the MSD boundary 



where Qb-i is the real-valued frequency of the solution near the boundary. Thus, Bb would have the form 
Bb = [h? /a)Q,b-i. Therefore, we can linearize the MSD boundary condition by treating the Bb term as a 
constant (which can change over the course of the simulation, similar to the nonlinearity of the NLSE). 
When inserted into the NLSE, the MSD boundary condition of Eq. (fT2|) yields 



VHb 



Im ( z -^^^ ) + - {Nb-i - Nb) 



*6 



where 



and therefore Dh 



Nb = s\^b\-Vb, 



Nb-^ = s I* 



6-11 



-Vb- 



(13) 



(14) 



Im 



0^fe^)+i(^^-i-^^^) 



This too is a nonlinear, non-constant term, and 



so must be treated as a constant in the same manner as the nonlinearity of the NLSE. 



3.4 Laplacian-Zero Boundary Condition 

The Laplacian-zero boundary condition is defined by 

and therefore Db — 0. In terms of the time-derivative of the NLSE, the LO boundary condition is given as 



5* 

'dt 



iis\^b\''-'Vb)^b 



making Bb — (/i^/a)(s|*6p — Vb). This condition is as easy to implement as the Dirichlet, and can be useful 
in many situations. 

To assist the stability analysis, a summary of the values of Bb and Db for all the mentioned boundary 
conditions are given in Table. [1] for future reference. Many other boundary conditions exist for simulating 
the NLSE, in which case the analysis shown in this paper can be adapted to the other boundary conditions. 



Table 1; Boundary condition terms for use with stability analysis. 



Boundary Condition 


Bb 


Db 


Dirichlet 







^{Vb-s\^b?) 




Laplacian-zero 


a 









MSD 


— Im 
a 


. *b-i . 




h^ 


Im 


I — r + - {Nb-i - 


-Nb) 





4 One-Dimensional Stability Analysis 

In the one-dimensional cases we analyze all four boundary conditions mentioned in Sec. [31 As stated, periodic 
boundary conditions yield a matrix where (in the linear case with s = and ^(r) = 0) the eigenvalues can 
be computed analytically. This allows the results obtained using the upper-bound methods (which we use 
with other boundary conditions) to be compared with the true eigenvalues giving an idea of how accurate 
they are. 

4.1 Second-Order Central Difference 

The second-order central difference is one dimension is given by 



V^*,; 



92^ 



dx^ 



^,+1 - 2^, + ^,_i 



and when implemented into the A matrix, forms a matrix which is tridiagonal (except for the two boundary 
condition rows). 



4.1.1 Periodic Boundary Conditions 

In order to obtain analytic expressions for the eigenvalues of A^ we start with the LSE case with no external 
potential and periodic boundary conditions. This yields the matrix 
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which, as per Def. [21 is a circulant matrix with c — {—2, 1, 0, ..., 0, 1}. Also, since A is a real-valued symmetric 
matrix, by Thm. [T] all eigenvalues are real and therefore the stability criteria of Eq. © can be used. By 
Thm. [31 the eigenvalues of A are given by 



exp 



2711 j 



N 



exp 



2TTij{N - 1) 



N 



j€{Q,...N~l}. 



The maximum value of \\j\ occurs either at j — N/2 if N is even, or j — [N ± l)/2 if N is odd. For N 
even- valued we have 

|A|max = -2 + exp [tti] -I- exp [tti] 



which yields 

For N odd-valued we have 

which yields 



= 4. 



\Mr 



-2-(-l)i/^ + (-l)^(-l)-i/^ 



|A|. 



-2-2cos(^) 



Taking TV — >■ oo, the maximum bound on the maximum absolute eigenvalue becomes 



|A|. 



<4. 



We therefore have an upper bound on the maximum absolute eigenvalue which, for even-valued iV, is guar- 
anteed to be one of the eigenvalues. The stability criteria of Eq. (jH]) is then formulated as 



k< 



VSh^ 



(15) 



In the general case where s 7^ and/or V{r) 7^ 0, the A matrix is no longer circulant (since the values 
of the nonlinearity or external potential vary over the diagonal of A). To get a bound on the maximum 
absolute eigenvalue, we make use of Thm. [51 The matrix A has N Gershgorin disks, since each diagonal 
entry of A can be unique, but each disk has the same radius (r = 2). Also, since the diagonal entries can in 
theory take on any value, all Gershgorin disk limits must be examined. This yields the stability bound 



k< 



V8 



h^ 



max{||L||oo,||i-4||oo} « 



(16) 



where we have defined the elements of L to be 



L, 



-is\^.\ 



V,) 



(17) 



where the index i spans over the entire grid. It is important to note that all values of L are 0{h?). Thus, 
for h ^ 1, and reasonable values of |^p and V, the linear bound of Eq. (ITSt should be very close to the true 
bound of the nonlinear problem. 

If we set L = in Eq. (fTB]). we recover the bound in Eq. (fT5|). This shows that (in this case at least), 
using the Gershgorin circle theorem yields the true bound on the eigenvalues of A. 



4.1.2 Dirichlet, MSD, and LO Boundary Conditions 

As shown in Sec. [31 Dirichlet, Laplacian-zero, and modulus-squared Dirichlet boundary conditions can all 
be viewed as single entries in the boundary value rows of the A matrix, denoted as B;,. As shown there, the 
values of Bh are real- valued and their values for each boundary condition were given in Table [H Using such 
a formulation, the A matrix becomes 
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In order to use the simple stability criteria of Eq. ©, we once again need to show that all eigenvalues of 

A are purely real. The A matrix is no longer symmetric, however it is easy to see that Bq and -Bat-i are 

eigenvalues of A, and the remaining eigenvalues of A are equivalent to the eigenvalues of the matrix A' 

defined as _ 

~ Li-2 1 
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Since A' is real- valued symmetric, we can use the stability bound of Eq. ([9]). 

We now need to find an upper bound on the absolute value of the eigenvalues of A' . We use the Gershgorin 
circle theorem to find all unique Gershgorin disks and take the limits of the disks to find the bounds on the 
absolute eigenvalues. Many of the Gershgorin disks are similar, differing only in the value of Li of the specific 
row. Therefore, each disk of different centers and radii has a subset of Li values relevant to it. Although in 
the current one-dimensional setting it is simple to define the subsets, in higher-dimensional settings, it can 
become burdensome to separate out each subset of L relevant to each Gershgorin disk of the same center 
and radius. Therefore, for practicality purposes, we define our bounds using all possible values of Li for 
each Gershgorin disk center and radius. This may make the resulting stability bound slightly higher than 
necessary in certain cases, but this is outweighed by the easc-of-use of the simplified bounds. The unique 
forms of the Gershgorin disks of A' are shown in Table [2] The resulting general stability bounds are 

Table 2: Unique forms of the Gershgorin disk centers (an) and radii (r^) for the A' matrix of the one- 
dimensional second-order central difference scheme. 
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ri = T,i^T Wzjl 
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V8 



max{||B||oo,||VLi,L. -GIloo} « 
where B are all boundary condition values (in this case Bq and Bat-i), and G is defined as 

G = {4,3,I,0}. 



(18) 



(19) 



In general, all possible values of G must be taken into consideration since there is no theoretical restriction 
on what values L can take. However, in certain specific circumstances, some of the values of G can be ignored 
(for example, when s < and V{r) > 0, only the largest magnitude value in G is needed). 



4.2 Fourth Order Central Difference 

The standard fourth order central difference scheme is given by 



V^*, 



^2^ 



dx^ 



-* 



i+2 



16*,+i - 30*, + 16*,_i - *,_2 



12 /i2 



(20) 



The stability analysis follows directly from the second-order case. The only major difference is that since the 
fourth order stencil is five points wide, the grid points near the boundary may need special consideration for 
the different boundary conditions. For our purposes here, we use the two-step high-order compact (2SH0C) 
version of the fourth-order scheme as described in Ref. |13| , in which case the near-boundary points can be 
formulated by combining the two steps of the 2SH0C scheme. 



4.2.1 Periodic Boundary Condition 

In the periodic case, no special attention is needed near the boundaries, and the A matrix in the LSE case 
with no external potential is 
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which is a circulant matrix, and its eigenvalues are therefore 
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The maximum absolute value once again occurs occurs at either j — N/2 if N is even, or j — (N ± l)/2 if 
N is odd. For N even-valued we have 



A 



15 4 1 1 -,\Af-2 I 4, ^^N-l 

^/^ = ^¥~3"12"12(~^) +3^"^^ 



16 



For N odd, we have 

A(^+i)/. = -'4-i f (-1)^/"^ + i-i)-'^'') i f (-1)^/^ 

which yields 
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'2/N 
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iN+l)/s 



15 4 / /7r\\ 1 / /27r 

-Y-3(2^°H]v))-12 2^°H^ 



As iV — > oo, |A| — ?► ^, which is the same bound as the iV-even case. Thus, the stability bound is given by 



k< 



3\ a/8/i^ 



(21) 



which we note is only a 25% reduction of the second-order bound of Eq. ([T5|) . In the general case where 
L y^ 0, the bound becomes 

k< ^V2 h^ 

max{||3L-16||oo,||3L+l||oo} a 

If s < and V{r) > 0, the first term in the denominator is the maximum of the two terms, and the resulting 
stability bound is equivalent to that found in Ref. [B] for using the RK4 scheme and fourth-order spatial 
differencing with the coupled NLSE. 
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4.2.2 Dirichlet, MSD, and LO Boundary Conditions 



As per Sec. [31 we formulate all three boundary conditions in terms of a Bf, entry in the A matrix. As 
discussed, an important issue is that we need to handle the grid points near the boundary due to the width 
of the scheme. A common way of dealing with the closest-interior points is to compute the Laplacian at 
those points using second-order differencing, however this can lead to the overall scheme becoming second- 
order. However, since we are using the 2SH0C version of the fourth-order differencing, we can derive the 
closest-interior points which, if the assumptions of the chosen boundary conditions hold, should maintain 
fourth-order accuracy. In one dimension, the 2SH0C scheme is defined as [13] 



1) A = -^ (*»+! - 2*, + *,_i) , 



V^*, 



I- 



-(A+i + A-i). 



(23) 
(24) 



In the first step, the second-order Laplacian is computed with the chosen boundary condition applied to 
it. Next, the result is used to compute the fourth-order Laplacian. As mentioned in Ref. [T^, this two- 
step scheme is equivalent to the standard wide-stencil of Eq. (I^Ul) for the interior points. We use the form 
of Eq. (ITU) for the boundary conditions on the Laplacian, and after combining the steps of Ec^. ([25)1 and 
Eq. (p4)) . we get the matrix 
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The Bb and Di, terms for each boundary condition are once again given in Table [TJ As in Sec. 14.1.21 the A 
matrix is not symmetric and has eigenvalues equal to B. The remaining eigenvalues are those of the matrix 
A' defined as 



A' 



^1~ T2 

4 
3 

_ J_ 

12 








is 



4 
3 


1 

12 








15 
6 


4 
3 


1 

12 





4 
3 


r 15 


4 
3 


1 

12 





1 
12 


4 
3 


Ln-a - 








1 
12 


4 
3 











1 

12 



15 



L 



15 
N-3 - -fi- 



_J_ 
12 

4 
3 



Lm- 



N-2 



29 
12 J 



which is once again real-symmetric so the bounds of Eq. ^ can be used. It is interesting to note that the 
values of Z?f, do not appear in any of the eigenvalues of A' . 

The unique forms (see the discussion in Sec. I4.1.2|) of the Gershgorin disk centers and radii of A' are 
shown in Table [31 The full stability bound is the same form as Eq. (|18l) , but with G defined as 



G=^x {64,63,46,12,-3,-4}. 



(25) 
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Table 3: Unique forms of the Gershgorin disk centers (an) and radii (r^) for the A' matrix of the one- 
dimensional fourth-order 2SH0C scheme. 



Am 


n = Y.^^iWij\ 


Li - 5/2 
Li - 5/2 
Li - 29/12 


11/4 
17/6 
17/12 



Once again, in the most general case, ah values of Eq. 
allowed time-step value. 



251) must be considered in finding the maximum 



5 Two-Dimensional Stability Analysis 

In higher dimensions, the A matrix is formed by unwrapping the solution into a one-dimensional vector and 
then formulating the scheme matrix accordingly. 

In Sections 14. 1.11 and 14. 2. II we noted that the stability bounds given using the Gershgorin circle theorem 
were equivalent to those obtained analytically for the linear case with periodic boundary conditions. We 
therefore justify relying exclusively on the Gershgorin theorem for higher dimensions, and focus on the 
stability bounds for Dirichlet, MSD, and Laplacian-zero boundary conditions (since the periodic boundary 
condition bounds will be a subset of the bounds computed for the other boundary conditions) . 

5.1 Second-Order Central Difference 

The second-order central difference scheme in two dimensions is given by 



^''^^,J = 



92* 



dx^ 






1 




1 




1 


-4 


1 




1 





1-,; 



(26) 



The corresponding A matrix has a tri-banded structure, with diagonal sub-sections corresponding to the 
boundary values. The form of the A matrix is shown in Fig. [5] (we do not show the values of the entries of 
the matrix due to space considerations, but they can be obtained through symbolic math codes). As in the 
one-dimensional case, all diagonal entries (which are the boundary value entries Bi,) of A are eigenvalues, 
and the remaining eigenvalues are real and equivalent to those of a matrix A' which is real-symmetric, thus 
allowing the use of the bounds in Eq. ©. The form of A' is also shown in Fig. [51 

The unique forms of the Gershgorin disk centers and radii for A' are shown in Table HI The stability 

Table 4: Gershgorin disk centers and radii for the A' matrix of the two-dimensional central-difference scheme. 



an 


n = Ei^, kijl 


L,-A 
L^-A 
L^-A 


2 
3 

4 



bounds are once again the same as in Eq. p^ with B being the set of all boundary values B^ and with G 
now defined as 

G = {0,1,2,6,7,8}. (27) 

In the linear case (s = 0) with no external potential and periodic, Dirichlet or Laplacian-zero boundary 
conditions, we get the linear stability bound 



8 a 



(28) 
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Figure 2: (Color online) Form of scheme matrix A and A' for second-order central differencing of the two- 
dimensional NLSE. The dots represent non-zero entries of the matrices. The matrices shown are for a 7 x 7 
grid. 



As before, since within the A matrix, all boundary, potential, and nonlinear terms are 0{h^) the simple 
bound of Eq. (|28p with slight adjustment can be used in many applications. 

5.2 Fourth-Order Central Difference 

The fourth-order central difference scheme in two dimensions is given by 



V"*, 



1 



12 h^ 







1 










-16 
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-16 


60 


-16 


1 






-16 










1 







*, 



The low-storage version of the 2SH0C equivalent scheme is defined as [T3] 



1) A. = ^ 





1 




1 


-4 


1 




1 





(29) 



(30) 



2) V^^.j « - 
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1 




1 


-12 


1 




1 





D 



«j 



1 



1 




1 




-4 




1 




1 



* 



«J- 



(31) 



The corresponding A matrix has a five-banded structure. The structure of the A matrix and its corresponding 
A' matrix are shown in Fig. [31 

The resulting unique forms of the Gershgorin disk centers and radii for A' arc shown in Table [5] The 
linearized stability bound is once again Eq. ([TS]) but with G defined as 

"^ X {128, 127,126, 110,109,92,24,9,8,-6,-7,-8}. (32) 



G 
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The linear bound (with s = and V{r) ~ 0) is then given by 

3%/8 }? 
32 a 



(33) 
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Figure 3: (Color online) Form of scheme matrix A and A' for the fourth-order 2SH0C scheme of the two- 
dimensional NLSE. The dots represent non-zero entries of the matrices. The matrices shown are for a 7 x 7 
grid. 



Table 5: Gershgorin disk centers (an) and radii (r^) for the A' matrix of the two-dimensional 2SH0C scheme. 



azi 


n^J2^^iWij\ 


Li -5 


11/2 


Li -5 


67/12 


Li -5 


17/3 


Li - 29/6 


17/6 


Li - 59/12 


25/6 


Li - 59/12 


17/4 
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which, as in the one-dimensional case, is 25% lower than the second-order linear bound given in Eq. 



6 Three-Dimensional Stability Analysis 

For the stability analysis in three dimensions, the same procedure utilized in the two-dimensional case of 
Sec. [5] is used. 

6.1 Second-Order Central Difference 

The second-order central difference scheme in three dimensions is given by 



V'*.,j,fe 



92^ 



1 



92^ 



i-j,k 



dy^ 



d^-q/ 



i,j,k 



dz^ 



ij.k 











1 











^ij'+l,fe + 





1 




*»,j,fc + 








1 


-6 


1 




1 






1 











(34) 



* 



i,j-l,k 



and the structure of the corresponding A and A' matrix are given in Fig. |4l The unique forms of the 
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'."A 
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Figure 4: (Color online) Form of scheme matrix A and A' for the second-order central difference scheme of 
the three-dimensional NLSE. The dots represent non-zero entries of the matrices. The matrices shown are 
for a 5 X 5 grid. 

Gcrshgorin disk centers (an) and radii (r^) for A' are shown in Table [6l The stability bounds of Eq. (jTSl) in 

Table 6: Gcrshgorin disk centers {an) and radii (r^) for the A' matrix of the three-dimensional central 
difference scheme. 



au 


n^Y.i^^a^3\ 




3 

4 
5 
6 



this case has G defined as 



G = {12,11,10,9,3,2,1,0}. 



(35) 
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In the linear case (s — 0) with no external potential and periodic, Dirichlet or Laplacian-zero boundary 
conditions, the linear stability bound becomes 



12 a 

6.2 Fourth-Order Central Difference 

The fourth-order central difference scheme in three dimensions is given by 
1 



V^* 



12 h^ 



[*-.+2,j.fe + **-2j,fc + 'J'..,j+2,fe + 'J'*,j-2./c + *»,j,/c+2 + *ij,fc-2 



-16 {■9^+l,J,k + ^^^-l,j.k + *.,j + l,fc + ^^^J-l,k + *^,j,/c+l + *ij,fc-l) + 90 *,j- fe] 

The single-storage version of the 2SH0C equivalent scheme is defined as [13] 
1) Aj,fc - 



1 



/ 












1 




\ 









* 



i,j-\-l,k 





1 




*»j\fc + 
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-6 


1 




1 






1 











* 



i,j — l,k 



(36) 



(37) 



(38) 



2) V^^.j, fc 



1 
12 



/ 
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Di.j+i^k + 
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1 






1 
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\ 




1 
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D^,oM + 
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-10 


1 




1 






1 











D 



i,j — l,k 



^2,J"+1,/C + 



1 




1 


*i,i,fc + 




1 






-12 




1 




1 


1 




1 




1 





*i,7-l,fc 



(39) 



and the structure of the corresponding A and A' matrix are given in Fig. [5l while the unique forms of the 
Gershgorin disk centers and radii for A' arc shown in Table [T] The stability bounds are the same as in 
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Figure 5: (Color online) Form of scheme matrix A and A' for the fourth-order 2SH0C scheme of the three- 
dimensional NLSE. The dots represent non-zero entries of the matrices. The matrices shown are for a 7 x 7 
grid. 
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Table 7: Gershgorin disk centers and radii for the A' matrix of the three-dimensional 2SH0C scheme. 



aii 


^* = E»^, l«ul 


Li - 15/2 


33/4 


Li - 15/2 


25/3 


Li - 15/2 


101/12 


Li - 15/2 


17/2 


Li - 22/3 


67/12 


Li - 22/3 


17/3 


Li - 29/4 


17/4 


Li - 89/12 


83/12 


Li - 89/12 


7 


Li - 89/12 


85/12 



Eq. p8|) . but with G now defined as 



G = -^ X {192, 191, 190, 189, 174, 173, 172, 156, 155, 138, 36, 21, 20, 6, 5, 4, -9, -10,-11, -12} . (40) 

In the linear case (s = 0) with no external potential and periodic, Dirichlet or Laplacian-zero boundary 
conditions, we get the linear stability bound 



, V8/i2 
16 a 



(41) 



which, as in the one- and two-dimensional cases, is simply 3/4 of the second-order bound given in Eq. ((36 



7 Conclusion and Summary of Results 

In this paper we have formulated linearized stability bounds for using second- and fourth-order spatial finite- 
differencing with fourth-order Runge-Kutta time-stepping for the multi-dimensional nonlinear Schrodinger 
equation (NLSE) with Dirichlet, modulus-squared Dirichlet, Laplacian-zero, and periodic boundary condi- 
tions. 

A summary of the stability results for easy reference is given presently. For the nonlinear Schrodinger 
equation defined as 

i ^— -H a V^* - y(r)* + s |*P^ = 0, 
at 

where a > and s are parameters of the system and ^(r) is an external potential, the numerical stability 
bounds on the time-step when using the fourth-order Runge-Kutta time-stepping scheme is as follows: 

In the linear case where s = and with no external potential (T^(r) = 0), utilizing periodic, Dirichlet, 
or Laplacian-zero boundary conditions, the stability bound on the time-step k when using the second-order 
central difference (CD) scheme in a d-dimensional setting is 



kcD < 



h^ 



dV2a 



(42) 



while that of using a fourth-order central difference scheme (with interior points computed in the two-step 
high-order compact (2SH0C) methodology of Ref. [TB]) is 



k2SHOC < 



4/ dV2a 



(43) 
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The linearized stability bounds for the general NLSE are 

V8 



k< 



max{||B||oo,||VL„L,-G||oo} 



(44) 



where B are the boundary points as defined by Table [8] (or in the periodic case is ignored), the elements of 
L is defined as 



a 



V, 



where the index i spans the entire grid, and G is a set of values defined in Table |9l determined by the 
dimension and method being used. 



Table 8: Values of B in Eq. (04 





Dirichlet (*& = 


= const) 


Laplacian-zero (V'^^f, = 


= 0) 


MSD (|*b|^= const) 


Bb 





^(.|*.P-H) 


— Im 

a 







Table 9: Values of G in Eq. g!]). 



Scheme -> 


CD 0(/i^) 


2SH0C 0{h'') 


ID 
2D 
3D 


{4,3,1,0} 

{8,7,6,2,1,0} 

{12,11,10,9,3,2,1,0} 


^x {64,63,46, 

12,-3,-4} 

— X {128,127,126,110,109, 

92,24,9,8,-6,-7,-8} 

-^ X {192,191,190,189,174, 

173,172,156,155,138,36,21, 
20,6,5,4,-9,-10,-11,-12} 



We have found through numerical testing that to ensure stability in all dimensions for typical problems, 
the bounds must be lowered by about 10%-20% (most likely due to nonlinear effects). Also, we note that 
the reduced linear results are often similar to the full linearized bounds and can therefore be used as a good 
estimate of the stability bound. 

Acknowledgments 

This research was supported by NSF-DMS-0806762 and the Computational Science Research Center (CSRC) 
at SDSU. We gratefully acknowledge insightful discussions with Peter Blomgren. 

References 

[1] L. Debnath. Nonlinear Partial Differential Equations for Scientists and Engineers. Birkhauser Boston, 
second edition, 2005. 



18 



[2] P.G. Kevrekidis, D.J. Frantzeskakis and R. Carrctero-Gonzalez. Emergent Nonlinear Phenomena in 
Bose-Einstein Condensates: Theory and Experiment, volume 45. Springer Series on Atomic, Optical, 
and Plasma Physics, 2008. 

[3] J.C. Butcher. A history of Runge-Kutta methods. Appl. Num. Math., 20, 3 (1996) 247-260. 

[4] J. C. Strikwerda. Finite Difference Schemes and Partial Differential Equations. SIAM, second edition, 
2004. 

[5] J. D. Lambert. Numerical Methods for Ordinary Differential Systems. John Wiley & Sons, 1991. 

[6] M. M. Cerimele, M. L. Chiofalo, F. Pistella, S. Succi and M. P. Tosi. Numerical solution of the Gross- 
Pitaevskii equation using an explicit finite-difference scheme: An application to trapped Bose-Einstein 
condensates. Phys. Rev. E., 62, 1 (2000) 1382-1389. 

[7] W. Dai. An unconditionally stable three-level explicit difference scheme for the Schrodinger equation 
with a variable coefficient. SIAM Jour. Num. Analy., 29, 1 (1992) 174-181. 

[8] M. S. Ismail. A fourth-order explicit schemes for the coupled nonlinear Schrodinger equation. Appli. 
Math. Comp., 196 (2008) 273-284. 

[9] R. Bronson. Schaum's Outline of Matrix Operations. McGraw-Hill, first edition, 1988. 

[10] G. H. Golub. Matrix Computations. Johns Hopkins University Press, third edition, 1996. 

[11] I. Kra and S. R. Simanca. On circulant matrices. Not. AMS, 59, 3 (2012) 368-377. 

[12] R. M. Caplan and R. Carretero-Gonzalez. A modulus-squared Dirichlet boundary condition for time- 
dependant complex partial differential equations and its application to the nonlinear Schrodinger equa- 
tion. Submitted to Appl. Math and Comp., (2012) arXiv:1110.0569. 

[13] R. M. Caplan and R. Carretero-Gonzalez. A two-step high order compact scheme for the Laplacian 
operator and its implementation in a fully explicit method for integrating the nonlinear Schrodinger 
equation. Submitted to Appl. Math and Comp., (2012) arXiv:1109.1027. 



19 



